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1. Introduction 

The PACS-CS (Parallel Array Computer System for Computational Sciences) project[jl], ^, 
|3|, ^ ^, ^, ^ |8]] aims at A^y^ = 2 + 1 lattice QCD calculations at the physical point to remove the 
most troublesome systematic errors associated with the chiral extrapolations. So far our simulation 
points cover from 67MeV to 3.5MeV for the degenerate up-down quark mass with the strange 
quark mass fixed around the physical value. The reduction of mud down to 10 MeV is achieved 
by the domain-decomposed Hybrid Monte Carlo (DDHMC) algorithm with the replay trick[^, p^ ]. 
For the simulation at mud = 3.5 MeV we incorporate some algorithmic improvements such as the 
mass preconditioning[ll , 12], the chronological inverter[13] and the deflation teqnique[p^] which 



make simulations stable and contribute to reduce the simulation cost. For the strange quark part we 



employ the UV-filtered Polynomial Hybrid Monte Carlo (UVPHMC) algorithm[115|]. 

In this report we present the simulation details and some eminent results for the hadron spec- 
trum. Chiral analyses on the pseudoscalar meson sector with the SU(2) and SU(3) chiral pertur- 
bation theories and calculation of the charm quark systems with the relativistic heavy quark action 



are given in separate reports! 16, 17 1. 



2. Simulation details 

We employ the C? (a) -improved Wison quark action with a nonperturbative improvement co- 



efficient Csw = 1.715[|18p and the Iwasaki gauge action[|19|]. All the simulations are carried out on 
a 32^ X 64 lattice at j8 = 1.90 corresponding to the lattice spacing of a = 0.09 fm. Table [T] sum- 
marizes our simulation parameters. We choose combinations of the hopping parameters (jCud, kTs) 
based on the previous CP-PACS/JLQCD results[|g, 0] except {K^d^K^) = (0.137785,0.13660) 
which is adjusted at the physical point with the use of the PACS-CS results in an early stage[^, ^. 
The physics results at (jCud, K"s) = (0.137785,0.13660) is presented in Ref. [Q]. 

The DDHMC algorithm is implemented for the up-down quark by domain-decomposing the 
full lattice with an 8^ block size as a preconditioner for HMC. The domain-decomposition factor- 
izes the up-down quark determinant into the UV and the IR parts geometrically. As a result we 
have the gauge force and the up-down quark force with the UV and the IR parts in the molecular 
dynamics evolution. The reduction of the simulation cost is achieved by applying the multiple time 
scale integrator[|^] to these three forces. We find that the relative magnitude of the force terms is 
given as follows: 



uvll • M-riRi 



16:4:1, (2.1) 



where Fq denotes the gauge part and Fuv.ir for the UV and the IR parts of the up-down quark. The 
associated step sizes 5tg, 5tuV) are chosen such that 

5tg I |Fg 1 1 « 5tuv I li^uv 1 1 ~ 5tir I |FiR 1 1 . (2.2) 

These step sizes are controlled by three integers No,Ni,N2 as Stq = r/NoNiNj, 5tuv = '^/N\N2, 
5xiR = t/N2 with T the trajectory length. We fixNo = Ni =4 in our all simulations. The value of 
is adjusted to make the simulation stable. The threshold for the replay trick is chosen to be dH = 2. 
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Table 1: Summary of simulation parameters. Quark masses are perturbatively renormalized in the MS 
scheme at the scale of fi = I /a. The replay trick is applied for the case of dH > 2. MD time is the number 
of trajectories multiplied by the trajectory length t. CPU time for unit r is measured on 256 nodes of the 
PACS-CS computer. 
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For the strange quark we employ the UVPHMC algorithm, where the domain-deconposition is not 
used. The polynomial order A^poiy for the UVPHMC algorithm is adjusted to keep high accep- 
tance rate for the global MetropoUs test at the end of each trajectory. Based on our observation 
of ||Fs|| II^irII for the strange quark force, we set St^ = Stir. Calclation of the IR force re- 
quires the inversion of the Wilson-Dirac operator on the full lattice, which is carried out by the 
SAP (Schwarz alternative procedure) preconditioned GCR algorithm. We use the SSOR precon- 
ditioned GCR algorithm for the UV part. These preconditionings are accelerated with the single 
precision arithmetic. We employ the stopping condition \Dx — b\/\b\ < 10~^ for the force calcula- 
tion and 10^'"^ for the Hamiltonian, which guarantees the reversibility of the molecular dynamics 
trajectories to high precision. The DDHMC algorithm for the up-down quark works efficiently for 
K-ud < 0.13770. 

As we reduce the up-down quark mass, the increasing fluctuations of the | \Fir\ \ make the sim- 
ulation unstable. To suppress the fluctuations of ||^ir||, we incorporate the mass preconditioning 
for the IR part (MPDDHMC), which spUts the IR force Fir into and Fir by introducing a new 
hopping parameter k^^ = Pi^iid with pi less than unity. In the MPDDHMC algorithm we need 
four integers {No,Ni,N2,N3) to controU the four step sizes 5tg, 5tuv, 5Tj'j^, Stjr. N2,Ni, and pi 
are adjusted to reduce the fluctuations of ||F/jj|| and ||Fir||. We choose Sz^ = for the strange 
quark force in the UVPHMC algorithm. 

For the run at JCud = 0.137785 further mass preconditioning is applied to the shifted IR force 
Fj^, which is divided into and Fj^ using an additional hopping parameter K^^ = p2 Kj^^j = piPi Kud 
with p2 less than unity. We refer to this algorithm as MP2DDHMC because of two-level of mass 
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preconditioning. In this case five step sizes 5tg, Stuv, 5Tj'^, 5fiR are controlled by five inte- 
gers {Nq ,Ni,N2,Nt,,N4). We adj ust tlie values of N2,Nt,, N4 and pi , p2 to keep stable the fluctuations 
of I I^irI U I^K I U l^iR W-Sr, is equal to 5^1'^. 

For the MPDDHMC and the MP2DDHMC algorithms the inversion of the Wilson-Dirac oper- 
ator on the full lattice is composed of three steps. Firstly, we prepare the initial solutions employing 
the chronological guess with the last 16 solutions. Secondly, we apply a nested BiCGStab solver 
consisting of the outer solver and the inner one. The latter with single precision arithmetic works 
as a preconditioner for the former operated with double precision. We employ a stringent stop- 
ping condition \Dx — b\/\b\ < 10^^^ for the outer solver and an automatic tolerance control ranging 
from 10"^ to 10"^ for the inner solver. Thirdly, the nested BiCGStab solver is replaced by the 
GCRO-DR (Generalized Conjugate Residual with implicit inner Orthogonalization and Deflated 
Restarting) algorithm, once the inner BiCGStab solver becomes stagnant during the inversion of 
the Wilson-Dirac operator. 

In Figs. H and ^ we show the dH and the force histories at JCud = 0.13727 and 0.13770 with 
the DDHMC algorithm and those at K^d = 0.13781 with the MPDDHMC algorithm. The time 
histories at JCud = 0.13727 are quite stable, whereas the K"ud = 0.13770 case shows the spike-like 
fluctuations of the IR force at a few % rate of trajectries. For the K"ud = 0.13781 run we observe 
that the MPDDHMC algorithm succeeds in reducing the fluctuations of the IR forces F/j^ and Fir. 




Figure 1: dH histories for {Kud,K,) =(0.13727,0.13640), (0.13770,0.13640) and (0.13781,0.13640) from 
the left. 



.1 .„..., k,., iiiJi. ii 



Figure 2: Force histories for (fCud, K,) =(0.13727,0.13640), (0.13770,0.13640) and (0.1 378 1,0. 13640) from 
the left. In the left and middle figures black, red and green lines denote Fq, Fuy and Fir, respectively, with 
the DDHMC algorithm. In the right figure black, red, green and blue lines are for Fq, F{jy, F/j^ and Fr, 
respectively, with the MPDDHMC algorithm. 



4 



2+1 flaavor lattice QCD simulation with 0{a)-improved Wilson quarks 



N. Ukita 



3. Hadron spectrum 

We measure hadron correlators at ever 10 trajectories for K^^ < 13770 and 20 trajectories for 
^ud > 13781. Light hadron masses are extracted from single exponential fits to the correlators 
with an exponentially smeared source and a local sink. In order to increase the statistics we take 
four source points with different time slices for K"ud > 0.13754. They are averaged on each configu- 
ration before the jackknife analysis. This reduces the statistical errors by typically 20-40% for the 
vector meson and the baryon masses and less than 20% for the pseudoscalar meson masses com- 
pared to a single source point. Statistical errors are estimated by the jackknife method. Figure || 
shows the binsize dependence of the error for the pion mass and the "t]v.«" meson mass. We ob- 
serve that the magnitude of the error reaches a plateau after 100-200 MD time. This feature seems 
almost independent of the quark mass. Since similar binsize dependences are found for other par- 
ticle types, we choose a binsize of 250 MD time at K"ud < 0.13770. At KTud = 0.13781 we employ a 
binsize of 1 10 MD time due to the lack of statistics. 
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Figure 3: Binsize dependence of the magnitude of error for rrijt (left) and m^q^^ (right) at K^d > 0.13754 
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Figure 4: Effective masses for the mesons (top) and the baryons (bottom) for (jCud, K^) =(0. 13727,0. 13640), 
(0. 1 3770,0. 1 3640) and (0. 1 378 1 ,0. 1 3640) from the left. 
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Figure 5: Light hadron spectrum extrapolated to the physical point (red circles) in comparison with the 
experimental values (black bars). 



Figure^shows the hadron effective masses at JCud = 0.13727,0.13770 and 0.13781. We ob- 
serve clear plateau for the mesons except the p meson at K"ud = 0.13781 and also good signal for the 
baryons thanks to a large volume. Especially, the Q. baryon has a stable signal and a weak up-down 
quark mass dependence for our simulation parameters. Taking advantage of this virtue we choose 
the Q. baryon as input to determine the lattice cutoff. Combined with the additional inputs of m^j 
and niK to determine the physical up-down and strange quark masses, we obtain a^^ = 2.176(31) 
GeV. In this procedure we employ the SU(2) ChPT analyses for the quark mass dependences of m^i, 
niK, fn and fx taking account of the finite size corrections evaluated at the one-loop level||8[ p^ ]. 
For niQi we assume the linear quark mass depenedences. With the use of this cutoff we find that 
the lightest pseudoscalar meson mass we have reached is about 160 MeV. To obtain the vector me- 
son masses and the baryon masses at the physical point we avoid the chiral analyses based on the 
heavy meson effective theory or the heavy baryon ChPT because of their poor convergences in the 
chiral expansions. We instead use linear chiral extraporations to the physical point. In Fig. |5| we 
compare the light hadron spectrum at the physical point with the experimental values. The largest 
discrepancy is at most 3%, albeit errors are still not small for the p meson, the nucleon and the A 
baryon. It should be also noted that our results contain possible 0((fl'AQCD)^) cutoff errors. 

We calculate the bare quark masses using the axial vector Ward-Takahashi identity (AWI) 
defined by am^^^ = hm^^oo (V4A™P(f)/'(0))/(2(P(f)P(0))) where A™p is the nonperturbatively 
O (a) -improved axial vector current p3|]. Employing the perturbative renormalization factors Za 



and Zp evaluated up to one-loop level 1 24, 25 1, we obtain 



mf^i^l = 2GeV) = 2.527(47)MeV, mf ^(^t = 2GeV) = 72.72(78)MeV. (3.1) 

The physical up-down quark mass is 30% smaller than our lightest one = 1 /a) = 3.5 MeV 

at (jfud, K^s) = (0.13781,0.13640). The results for the pseudoscalar meson decay constants are given 
by 

A = 134.0(4.2)MeV, = 159.4(3. l)MeV, /^/^ = 1.189(20) (3.2) 

at the physical point with the perturbative Za- They are consistent with the experimental values 
within the errors. Our concern about the values for the quark masses and the pseudoscalar meson 
decay constants is the use of the perturbative renormalization factors which might cause sizable 
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systematic errors. We are now calculating the nonperturbative Za and Zp with the Schrodinger 
functional scheme. 
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